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Moment  Preserving  Adaptive  Particle  Weights  using  Octree 
Velocity  Distributions  for  PIC  Simulations 

Robert  Scott  Martin*  and  Jean-Luc  Cambier1 
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^ Spacecraft  Propulsion  Branch, 

Air  Force  Research  Laboratory, 

Edwards  AFB,  CA  93524  USA 

Abstract.  The  ratio  of  computational  to  physical  particles  is  of  primary  concern  to  statistical  particle  based  simulations  such 
as  DSMC  and  PIC.  An  adaptive  computational  particle  weight  algorithm  is  presented  that  conserves  mass,  momentum,  and 
energy.  This  algorithm  is  then  enhanced  with  an  octree  adaptive  mesh  in  velocity  space  to  mitigate  artificial  thermalization. 
The  new  octree  merge  is  compared  to  a  merge  that  randomly  selects  merge  partners  for  a  bi-Maxwellian  velocity  distribution. 
Results  for  crossing  beams  in  a  fixed  potential  well  along  with  an  electrostatic  PIC  version  with  and  without  MCC  collisions 
based  ionizing  breakdown  show  the  advantages  of  the  merge  algorithm  to  both  fixed  particle  weights  and  randomly  selected 
merge  partners. 
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INTRODUCTION 

The  ratio  of  computational  to  physical  particles  is  a  key  factor  in  determining  the  statistical  scatter  and  accuracy  in 
particle-based  simulation.  This  is  particularly  true  for  problems  characterized  by  wide  ranges  of  number  density  such 
as  those  found  in  spacecraft  electric  propulsion  plumes  and  ionizing  discharges,  where  populations  of  electrons  and 
excited  states  can  grow  exponentially.  A  particle  management  method  must  then  be  devised  which  balances  statistical 
accuracy  requirements  with  prevention  of  runaway  computational  costs. 

The  standard  approach  of  merging  of  particles[l]  using  pair-wise  coalescence  (2:1  ratio),  cannot  guarantee  simul¬ 
taneous  conservation  of  mass,  momentum  and  energy.  The  pair-wise  coalescence  conserves  mass  by  simply  summing 
the  computational  weight  of  the  individual  constituents.  One  computational  particle  only  has  three  degrees  of  freedom 
for  velocity.  The  resulting  particle  can  therefore  either  match  the  average  momentum  or  the  kinetic  energy  of  the  orig¬ 
inal  pair,  but  never  both  simultaneously.  The  growth  of  this  error  is  limited  by  selecting  particles  that  are  near  each 
other  in  velocity  space,  but  is  a  fundamental  consequence  of  the  reduction  of  degrees  of  freedom. 

As  a  result  of  this  inherent  limitation,  various  sophisticated  models  have  been  designed  to  mitigate  the  error 
in  momentum  or  energy.  Hewett  used  computational  particles  with  internal  energy  in  which  the  error  could  be 
accumulated[2].  Assous  and  then  Welch  designed  methods  of  merging  values  to  grid  nodes  and  redistributing  the 
moments  to  particles [3,  4].  Though  exact  energy  and  momentum  conservation  is  possible  with  these  methods,  they  are 
considerably  more  complicated  than  the  original  naive  2:1  merge. 


MOMENT  PRESERVING  MERGE 

Having  identified  the  lack  of  sufficient  degrees  of  freedom  in  the  merge  resultant  particle  as  the  source  of  the  error, 
Cambier  devised  a  simple  method[5]  which  relies  on  the  generation  of  a  pair  of  particles.  The  particle  pair  provides  the 
required  freedom  to  conserve  all  moments  up  to  2nd  order  exactly.  The  pair  of  resultant  particles  have  the  same  mean 
momentum  as  the  originals,  but  they  also  have  equal  and  opposite  components  of  velocity  perpendicular  to  the  mean 
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momentum  such  that  energy  is  conserved.  The  moments  to  be  conserved  are  shown  in  Equation  1 .  In  these  equation, 
the  ’k’  original  particles,  denoted  with  superscript  ’(/?)’,  and  the  bar  denotes  an  average  quantity  used  in  the  Merge. 
The  subscripts  denote  the  vector  direction  whether  index,  i,  or  directions  parallel  or  perpendicular.  Assuming  the  pair 
of  merge  result  particles  are  equal  weight,  the  formulas  for  the  merge  can  be  seen  in  Equation  2  where  superscript 
’(a /by  refers  to  one  or  the  other  of  the  merge  result  particles.  It  can  be  shown  that  the  moment  sums  are  equivalent 
between  the  k  original  and  pair  of  resultant  particles.  Thus,  pair-wise  reduction  is  obtained  through  an  equivalent  ratio 
of  4:2.  These  formulas  also  work  for  arbitrary  ratios  (n+2):2  with  the  same  conservation  properties. 


m  =  (1) 
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Vi~  I 

E*m^(v \p)-vl)2 

V ■  =  - - - — 

It  is  important  to  note  that  in  this  example,  some  ambiguity  in  the  sign  of  perpendicular  velocity  component  remains. 
For  the  purposes  of  this  work,  each  cardinal  direction  is  randomly  assigned  a  sign  for  one  of  the  particles  and  the 
second  particle  uses  the  opposite  sign  in  each  direction.  As  a  result  of  this  choice  of  signs,  though  the  component 
of  temperature  is  conserved  in  each  direction  independently,  though  the  choice  is  not  unique  and  particles  may 
scatter  perpendicular  to  a  dominant  axis  of  the  distribution.  A  better  choice  would  be  to  determine  the  sign  of  mixed 
second  order  moments  and  use  that  to  assign  the  signs  for  the  particle.  If  for  example,  most  of  the  particles  lay  on 
the  (—  v,  —  y,  —  z)  — »  (+v, +y, +z)  line,  the  sign  of  the  mixed  moments  would  all  be  positive  and  all  of  the  particle 
signs  should  match.  However,  this  approach  is  left  to  future  consideration  along  with  higher  moment  conservation  in 
extension  and  conclusion  section. 

Though  this  merge  exactly  conserves  the  moments  through  energy,  it  can  result  in  artificial  thermalization  of  the 
velocity  distribution.  This  is  easily  seen  by  considering  the  result  of  merging  particles  from  two  mono-energetic  beams 
of  opposite  direction.  Though  the  energy  of  the  resulting  velocity  distribution  would  be  unchanged,  the  merge  results 
in  a  range  of  particles  scattered  throughout  the  gap  between  and  outside  the  original  beams.  For  example,  if  two  equal 
mass  particles  were  selected  from  the  +v  side  beam  and  a  third  of  equal  mass  from  the  —  v  side,  the  mean  velocity 
would  be  at  +v/3  and  the  resulting  particle  velocities  would  be  symmetric  v/3  =b  y/24/27  v  and  therefore  clearly  not 
back  into  the  original  +v  and  — v  beams. 
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OCTREE  ADAPTIVE  MESH 

The  present  work  extends  this  exact  moment-preserving  merge  through  an  octree-based  adaptive  mesh  in  velocity 
space  to  ensures  that  merging  partners  are  relatively  close  in  phase  space.  This  mitigates  artificial  thermalization  due 
to  merging  of  particles  with  large  opposite  velocities  such  as  those  found  in  the  beam-beam  interaction  example. 
Rather  than  randomly  selecting  merge  partners,  the  particles  are  first  sorted  into  octants  of  velocity  space.  Each  of 
these  octants  is  then  recursively  subdivided  and  the  particles  are  sorted  into  the  sub-cells.  The  recursion  was  originally 
terminated  once  the  number  of  particles  in  a  given  sub-cell  is  below  a  specified  threshold.  The  conservative  merge  is 
then  applied  to  the  particles  within  one  velocity  space  sub-cell. 

To  test  the  modified  octree  based  conservative  merge,  a  bi-Maxwellian  velocity  distribution  was  used.  Figure  1 
contrasts  results  from  a  randomly  selected  4:2  merge  and  the  octree  based  (3-1 1):2  based  merge  where  (3-11)  refers 
to  the  range  of  computational  particles  per  sub-cell  to  be  merged.  In  the  example  shown,  the  random  merge  resulted 
in  a  reduction  of  10,000  — )>  5,000  computational  particles  while  the  octree  merge  starts  with  the  original  particles  and 
performs  a  reduction  of  10,000  -^4,727. 

The  threshold  was  set  to  1 1 -particles  per  sub-cell  because  it  was  found  to  result  in  a  net  merge  ratio  of  approximately 
2:1.  This  is  lower  than  5.5:1  because,  in  each  of  the  eight  children  cells,  the  number  of  particles  ranges  between  0- 
11  rather  than  being  exactly  11  particles.  For  example,  cells  with  fewer  than  three  particles  cannot  be  conservatively 
merged  at  all.  The  effective  merge  was  found  heuristically  to  result  in  approximately  2: 1  reduction  had  1 1  as  the  upper 


Random  Merge 


Octree  Merge 


-1500  -1000  -500  0  500  1000  1500  -1500  -1000  -500  0  500  1000  1500 

Vx  (m/s)  Vx  (m/s) 


FIGURE  1.  Comparison  of  randomly  selected  4:2  merge  (left)  and  octree  based  (3-ll):2  merge  (right).  The  scatter  plots  (top) 
show  particle  positions  in  velocity  space  and  the  line  plots  (bottom)  compare  the  original  and  merged  velocity  distributions 
integrated  over  y-  and  z-velocity.  In  the  scatter  plots,  the  original  particles  are  marked  by  the  black  dots,  the  left  and  right  halves  of 
the  original  bi-Maxwellian  are  circled  with  blue  and  red  respectively,  and  the  output  merged  particles  are  marked  in  violet. 


limit.  This  ratio  may  be  problem  dependent  and  should  therefore  be  considered  a  tunable  parameter  that  quantifies  the 
aggressiveness  of  the  merging  procedure  rather  than  as  a  direct  prescription  for  a  specific  merge  ratio. 

For  fully  adaptive  particle  weights,  an  analogous  particle  split  method  is  necessary  to  re-populate  depleted  velocity 
distributions  that  result  from  the  particle  merging.  Without  this  split,  particles  merged  when  the  spatial  distribution  is 
compact  result  in  extremely  noisy  values  as  the  distribution  expands.  The  split  is  achieved  by  performing  operations 
using  the  conservative  moment  sums  similar  to  those  of  the  merge.  Instead  of  consuming  the  entire  original  particles 
to  form  the  merge  results,  only  a  fraction  of  each  original  particle  is  consumed  to  form  the  new  pair.  The  effective 
split  has  a  ratio  of  (n+2):(n+2+2).  This  approach  only  offers  a  weak  growth  in  the  number  of  computational  particles. 
Instead,  subsets  of  3-4  particles  are  selected  from  the  particles  residing  in  the  particular  velocity  sub-cell  in  which 
the  split  operation  is  being  performed.  For  each  of  these  subsets,  a  3:5  or  4:6  split  is  performed  yielding  a  higher 
computational  particle  growth  rate.  This  split  technique  is  used  in  the  results  presented. 

In  later  studies  with  DSMC  collisions,  low  density  regions  were  populated  by  collision  events  that  scattered  velocity 
out  of  the  original  distributions.  Repeated  merging  as  particles  scattered  out  of  the  core  of  the  distribution  resulted  in 
a  depletion  of  computational  particle  number  within  the  core  in  favor  a  few  heavy  particles.  To  compensate  for  this 
trend,  the  cell  merge  limit  was  modified  such  that  cells  are  only  subdivided  when  the  total  weight  of  particles  within 
the  cell  exceed  a  target  weight.  This  target  weight  can  then  be  modified  for  each  cell  so  that  the  merge  procedure 
attempts  to  achieve  a  uniform  number  of  computational  particles  per  cell  and  therefore  computational  work. 
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FIGURE  2.  Position  vs  Time  results  of  particle  number  density  for  Analytical,  Fixed  Weight,  Random  4:2  Merge,  and  Octree 
based  Merge. 


POTENTIAL  WELL  TEST  CASES 

The  combined  fully- adaptive  particle  weighting  scheme  can  now  be  applied  to  several  test-problems.  The  first  is 
collisionless  thermal  beams  in  a  potential  well.  The  initial  conditions  consist  of  a  bi-Max wellian  distribution  of 
particles  at  the  center  on  a  fixed  parabolic  potential  well.  The  two  halves  of  the  bi-Max  wellian  initially  separate  in 
opposite  directions.  The  acceleration  due  to  the  electric  field  then  decelerates  the  particles  and  they  fall  back  across 
the  bottom  of  the  well.  They  then  proceed  to  the  opposite  side  of  the  well  and  continue  oscillating  indefinitely. 

This  setup  is  a  particularly  useful  surrogate  for  the  expansion  of  a  plasma  plume  because  the  particle  density  falls  off 
by  several  orders  of  magnitude  between  the  center  and  edges  of  the  well.  It  is  also  useful  because  an  exact  analytical 
solution  can  be  derived  for  the  evolution  of  the  density  profile  as  follows. 

In  this  simple  invariant  potential  well  setup,  every  particle  in  the  initial  distribution  follows  a  simple  oscillatory 
trajectory  such  that  the  sum  of  kinetic  and  potential  energies  is  constant.  It  can  then  be  shown  that  the  position  and 
velocity  for  each  particle  follows  the  sinusoidal  trajectory  based  on  the  initial  velocity  and  position. 

The  trajectories  can  then  be  solved  for  the  initial  velocities  for  particles  between  the  ’inner’  and  ’outer’  edges  of 
a  spatial  cell  at  any  time  of  interest.  These  inner  and  outer  velocities  can  then  be  used  as  the  lower  and  upper  limits 
of  integration  on  the  initial  bi-Maxwellian  velocity  distribution  to  determine  the  number  of  particles  within  a  cell  at  a 
given  time. 

For  the  particle  simulations,  an  implicit  Crank-Nicolson  particle  push  using  the  analytically  defined  parabolic 
potential  was  used.  This  push  was  selected  due  to  its  exact  marginal  stability  and  energy  conservation  properties 
for  simple  harmonic  oscillation.  This  allows  the  simulation  to  run  for  several  beam  oscillations  without  numerical 
heating  or  cooling.  The  timesteps  were  selected  based  on  a  CFL  criteria  such  that  the  fastest  particles  in  the  original 
distribution  only  cross  a  single  computational  cell  in  a  timestep. 

Figure  2  shows  number  density  results  for  the  Analytical,  Fixed  Particle  Weight,  Randomly  Selected  4:2  Merge, 
and  Octree  Merge.  Clearly  the  results  for  Octree  Merged  and  Fixed  particle  weights  are  quite  similar  and  match  the 
analytical  results  across  the  4-orders  of  magnitude  in  particle  density  depicted. 

This  problem  explicitly  highlights  the  disadvantage  of  the  random  merge  partner  selection  in  thermalizing  the 
distribution.  The  thermalization  occurs  in  the  first  timestep  when  the  particles  all  reside  in  the  central  cell  with  opposite 


FIGURE  3.  Comparison  of  number  of  computational  particles  throughout  the  simulations  for  Fixed  Weight,  Random  4:2  Merge, 
and  Octree  based  Merge. 


directions.  The  random  selection  mixes  particles  from  opposite  sides  of  the  initial  bi-Maxwellian  distribution  and  fills 
in  the  center.  Even  though  the  figure  looks  dramatically  different,  all  of  the  moments  are  conserved  in  each  cell.  For 
every  particle  that  lost  momentum  and  energy  to  fill  in  the  center  figure,  others  gained  momentum  and  energy  which 
results  in  a  wider  fringe. 

Figure  3  shows  the  number  of  computational  particles  for  the  different  particle  management  techniques.  All  of  the 
cases  started  with  the  same  6000  computational  particles.  The  figure  shows  that  the  initial  merge  between  the  random 
and  octree  techniques  are  of  similar  magnitude  despite  the  dramatic  difference  in  the  density  plots.  The  random  merge 
technique  also  does  not  include  a  split  and  so  the  number  of  particles  monotonically  decreases  in  steps  with  the  bounce 
frequency. 


SPARK  GAP  TEST  CASE 

A  second  test  of  the  merging  algorithm  involves  electron  avalanche  in  a  spark  gap  gas  breakdown.  The  simulation  is 
initialized  with  a  6kV  potential  drop  across  a  5mm  gap.  In  this  example,  the  particles  tracked  are  the  electrons  in  a 
partially  ionized  background  Argon  gas  which  is  assumed  stationary  for  the  duration  of  the  simulation. 

The  1-D  code  was  modified  to  include  the  effect  of  the  charge  distribution  on  the  electrostatic  potential.  The  charge 
density  p,  was  assumed  constant  within  the  cell  and  was  simply  the  difference  of  electron  and  ion  densities  multiplied 
by  the  fundamental  electron  charge.  Control  points  for  a  quadratic  b-spline  were  then  solved  for  the  potential  from 
the  charge  density  such  that  V2T>  =  — p.  The  boundary  conditions  of  ±3kV  were  also  applied  to  the  solution  of  the 
control  points.  Electrons  that  arrived  at  the  +3kV  boundary  were  re-injected  at  the  -3kV  boundary  with  a  small  thermal 
velocity  to  maintain  charge  neutrality.  The  injected  particles  were  injected  randomly  in  the  inside  half  of  the  boundary 
cell  so  that  self-forcing  due  to  the  OD-representation  of  the  charge  density  does  not  cause  the  particles  to  be  pushed 
back  out  of  the  domain.  Though  this  is  only  a  rough  representation  of  an  electrostatic  PIC  model,  it  is  sufficient  for 
testing  the  octree  merge  algorithms. 

In  addition  to  the  modifications  to  the  potential  for  the  spark  gap  case,  the  push  was  modified  to  an  analytical  form 
rather  than  using  Crank-Nicolson.  Because  the  potential  is  represented  by  a  piecewise  quadratic  b-spline,  the  electric 
field  is  a  piecewise  linear  continuous  function.  The  exact  solution  for  a  particle  undergoing  an  acceleration  that  depends 
linearly  on  position,  a  =  %-E  =  ^  [Eq  +  ||  |Xo (x  —  xo)],  can  easily  be  found  in  ID.  It  has  the  form  shown  in  Equation  3 
with  constants  defined  in  Equation  4.  The  equation  assumes  a  piecewise-linear  electric  field  and  is  therefore  only  valid 
within  one  computational  cell. 
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FIGURE  4.  Position  vs  Time  results  of  electron  density  and  computational  particles  per  cell  in  spark  gap.  Electrons  are  pushed 
from  the  cathode  on  the  left  to  anode  on  the  right  of  the  figures  due  to  the  6kV  applied  potential.  The  left  figure  shows  the  original 
electron  density  with  fixed  particle  weight  and  the  right  figure  shows  the  electron  density  for  the  octree  merge  algorithm. 
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To  account  for  the  change  of  slope  at  cell  boundaries,  the  cell  boundary  crossing  times  are  found  and  particles 
are  pushed  only  within  one  cell  at  a  time.  The  time  at  which  a  particle  crosses  xo  or  xo  +  Ax  can  then  be  solved  by 
substituting  x(t)  =  xo  or  x0  +  Ax  and  a  =  into  the  formula  for  x(t)  and  solving  the  quadratic  equation  for  a.  The 
crossing  times  for  the  cell  edges  are  then  simply  t  =  VXln(a).  Note  that  each  side  has  two  possible  a-s  as  roots  of 
the  quadratic  equation.  The  set  of  four  crossing  times  may  then  each  be  positive,  negative,  or  imaginary.  Each  particle 
is  then  advanced  by  the  smallest  real  positive  time  or  full  timestep  if  shorter  using  Equation  3  repeatedly  until  the 
full  timestep  is  complete.  The  full  timestep  is  based  on  a  CFL  condition  assuming  an  electron  accelerated  by  a  6kV 
potential  cannot  cross  more  than  one  cell  per  timestep. 

In  the  test  case,  the  background  was  initialized  at  0.1%  singly  charged  ionization  in  a  le22/m3  background  of  Argon. 
The  evolution  of  the  electron  density  is  shown  in  Figure  4.  An  interesting  feature  is  the  oscillation  of  trapped  electrons 
that  are  pulled  back  from  the  initial  acceleration.  This  feature  is  a  result  of  the  potential  exceeding  the  anode  boundary 
condition  in  the  positive  column  followed  up  to  the  anode  sheath.  The  slowest  electrons  that  have  insufficient  kinetic 
energy  to  overcome  the  sheath  voltage  fall  back  towards  the  cathode  and  oscillate  indefinitely  because  the  plasma  is 
collisionless.  This  example  is  quite  similar  to  the  crossing  beams  in  the  potential  well  of  the  previous  example  and 
again  demonstrate  the  importance  of  the  octree  merge’s  minimal  thermalization. 

In  addition  to  the  modified  particle  push,  a  simple  inelastic  MCC  collision  process  was  added  between  the  electrons 
and  background  gas.  The  energy  dependent  cross  section  used  was  based  off  the  data  in  Reference [6]  using  the  form 
given  in  Equation  5.  In  the  equation,  C  is  the  slope  of  cross-section  vs  energy  above  the  ionization  threshold  energy, 
I.  The  maximum  cross  section,  omax  is  the  peak  cross  section  prior  to  a  further  reduction  with  increasing  collision 
energy.  For  Argon,  values  of  C  =  2e-17cm2/ev,  I  =  15.8ev,  and  omax  =  3.7e-16cm2  from  Reference  [6]  were  also 
used  to  approximate  the  electron  impact  ionization  cross  section.  For  this  sample  problem,  all  other  collision  types 
including  Coulomb  and  recombination  collisions  were  ignored. 
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FIGURE  5.  Position  vs  Time  diagram  of  electron  density  and  computational  particles  per  cell  in  an  ionizing  spark  gap  using 
MCC  collision  model.  Electrons  are  pushed  from  the  cathode  on  the  left  to  anode  on  the  right  of  the  figures  due  to  the  6kV  applied 
potential.  The  left  two  figure  shows  the  original  electron  density  and  computational  particles  per  cell  with  fixed  particle  weights 
and  the  right  two  figure  shows  the  same  results  using  the  octree  merge  algorithm. 


a(e)  =  rnin(C(e-I),omax)  (5) 

Based  on  the  collision  frequency  from  the  cross  section,  particles  are  selected  to  either  undergo  collision  with  the 
background  gas  or  not.  If  a  particle  is  selected  for  collision,  it’s  kinetic  energy  is  reduced  by  the  ionization  potential,  /, 
and  a  new  electron  of  equal  computational  weight  is  created  with  zero  velocity  at  the  same  location.  The  background 
neutral  density  is  also  decreased  by  the  same  weight  so  that  the  ion  density  can  be  determined  as  ni  =  no  —  nn. 

With  the  addition  of  MCC  ionization  collisions,  the  growth  of  the  number  of  computational  electrons  is  exponential 
from  a  small  initial  ionization  fraction.  This  is  exactly  the  type  of  situation  for  which  the  merge  algorithm  was  designed. 
Without  the  ability  to  adaptively  adjust  particle  weights  throughout  the  simulation,  the  computational  cost  rapidly 
exceeds  the  resources  available  even  if  only  a  few  and  therefore  statistically  inaccurate  seed  electrons  are  initialized. 

Finally,  Figure  6  shows  the  total  number  of  particles  for  the  fixed  weight  and  octree  merged  MCC  ionizing 
breakdown  case.  The  particle  weight  adaptation  is  clearly  compensating  for  what  would  otherwise  be  an  exponential 
growth  in  computational  particles.  Though  the  density  results  are  noisier,  this  should  be  expected  for  few  degrees  of 
freedom  in  the  system.  The  key  feature  is  that  the  computational  cost  is  now  decoupled  from  particle  weight  using  the 
octree  based  particle  merge. 


FUTURE  EXTENSION  AND  CONCLUDING  REMARKS 

The  issue  of  mixed  second  order  moments  mentioned  previously  ties  into  the  concept  of  extending  the  conservative 
merge  to  higher  moments  of  the  velocity  distribution.  Using  index  notation,  it  is  clear  that  current  formulation 
conserves  one  0th  order  moment  for  the  mass,  m,  three  1st  order  moments  for  velocity  vj,  and  three  diagonal  2nd  order 
moments  V;V;.  In  the  typical  notations  of  the  moments  in  statistical  mechanics,  it  is  useful  to  define  these  moments 
instead  in  terms  of  the  peculiar  velocity  of  the  particles,  Vi  =  v*  —  vj.  The  diagonal  second  order  moment  is  then  simply 
proportional  to  the  diagonal  terms  of  the  pressure  tensor,  P[j  —  nmViV )  as  defined  in  Reference  [7].  Though  the  net 


FIGURE  6.  Comparison  of  number  of  computational  particles  throughout  the  ionizing  breakdown  for  Fixed  Weight  and  Octree 
Merge. 


thermal  energy,  V/V*  is  conserved,  each  component  of  ViVi  =  VXVX  +  VyVy  +  VZVZ  is  also  conserved  independently.  This 
allows  for  anisotropic  pressure/temperature,  but  makes  no  guarantee  that  shear  stress  is  conserved. 

Adding  a  third  product  particle  to  the  merge  procedure  with  the  constraint  that  all  three  product  particle  masses  are 
equal  provides  sufficient  degrees  of  freedom  to  satisfy  conservation  of  the  full  second  order  moment,  V[Vj.  Though  it 
would  appear  that  there  are  6-additional  moments  for  the  full  second  order  tensor,  the  moment  integral  is  commutative 
for  ViVj  =  VjVi  and  therefore  the  three  additional  degrees  of  freedom  are  sufficient.  However,  these  additional  degrees 
of  freedom  no  longer  allow  for  direct  solution  and  require  a  non-linear  solve  for  the  merge  result  velocities. 

Similarly,  adding  a  fourth  particle  adds  the  ability  to  conserve  VfVjVj.  This  is  equivalent  to  conserving  the  heat  flux 
vector,  and  is  sufficient  to  conserve  the  13-moments  necessary  to  the  Navier-Stokes  solution  to  the  Chapmann-Enskog 
expansion  without  the  restriction  of  isotropy.  Work  has  begun  to  solve  for  the  set  of  four  particle  velocities  that  result  in 
the  correct  moments  using  a  Newton-Raphson  style  non-linear  solver.  The  solution  requires  stabilization  due  to  saddle 
points  in  the  12-dimensional  solution  space.  Further  work  to  extend  this  method  to  solve  on  the  energy  constrained 
manifold  is  hoped  to  alleviate  these  issues.  Developing  a  merge  routine  to  this  level  would  be  potentially  useful  in 
developing  hybrid  particle/continuum  codes. 

Along  with  these  velocity  moments,  spatial  moments  can  be  similarly  conserved  and  are  necessary  for  conserving 
electrostatic  potential  energy  during  the  merge  as  developed  in  the  original  development  of  the  conservative  merge 
algorithm  [5].  In  addition,  future  work  will  include  consideration  of  conservation  of  magnetic  moments  and  internal 
energy. 

These  merging  and  splitting  algorithms  show  tremendous  promise  for  future  particle  codes  that  involve  strong 
particle  growth  and/or  large  dynamic  range  in  particle  density.  They  also  provide  the  groundwork  for  modified  collision 
routines  that  are  not  constrained  to  constant  numbers  of  computational  particles  as  discussed  in  the  accompanying 
paper,  Reference  [8]. 
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©  Potential  Function  of  e~  and  Ar+ 


Unmerged 


Position  (mm) 


15  16  17  18 
Electron  Density,  Iog10(  ng  ) 
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Gas  Breakdown 


o  Merge  Needed  w/  Exponential  #  Growth 
©  Examples... 

Chain  Branching:  H2  +  M  — >  2H  +  M 
Ionization:  Ar°  +  e~  Ar+  +  e~  +  e~ 

©  Ionizing  Breakdown  in  6kV  Potential 
©  Electrons  Flow  from  Cathode  — »  Anode 
©  Inelastic  MCC  with  Background 
©  Potential  Function  of  e~  and  Ar+ 

©  Merge  Retains  Features  and  Magnitude 
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o  Merge  Needed  w/  Exponential  #  Growth 
©  Examples... 

Chain  Branching:  H2  +  M  — >  2H  +  M 
Ionization:  Ar°  +  e~  Ar+  +  e~  +  e~ 

©  Ionizing  Breakdown  in  6kV  Potential 
©  Electrons  Flow  from  Cathode  — »  Anode 
©  Inelastic  MCC  with  Background 
©  Potential  Function  of  e~  and  Ar+ 

©  Merge  Retains  Features  and  Magnitude 


While  Controlling  Computational  Cost 
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Higher  Order  Moment  Conservation 


0 


o  Moments  defined  as  Integrals  of  VDF:  Q  =  f  Qnfdvi 
o  Discrete  Version:  nf  w^S(v^)  such  that  Q  =  w^Q/  J2 
o  Merged  Particles  have  4  DOF  each:  w,  vx,  vz 
o  Number  of  Moments  Conserved  from  Number  of  DOF 


Cartesian  Moments 


Moment 

Order 
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Higher  Order  Moment  Conservation 


o  Moments  defined  as  Integrals  of  VDF:  Q  =  f  Qnfdvi 
o  Discrete  Version:  nf  w^S(v^)  such  that  Q  =  Y  w^Q/  Y 
o  Merged  Particles  have  4  DOF  each:  w,  vx,  vy,  vz 
©  Number  of  Moments  Conserved  from  Number  of  DOF 


Moment 

Order 

Mass 

Mass  Flux 

0th 

Ist 

Y  =  W 

Y  =W'Vj 

1  Particle  -  Mass  &  Momentum 
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Higher  Order  Moment  Conservation 


o  Moments  defined  as  Integrals  of  VDF:  Q  =  f  Qnfdvi 
o  Discrete  Version:  nf  w^S(v^)  such  that  Q  =  Y1  w^Q/  ^ 
o  Merged  Particles  have  4  DOF  each:  w,  vx,  vy,  vz 
o  Number  of  Moments  Conserved  from  Number  of  DOF 


Moment 

Order 

Mass 

Mass  Flux 

Momentum  Flux 

0th 

1st 

2nd 

Y  =  w 

=  w-Vi 

=w-m 

2  Particles  -  Mass,  Momentum,  and  Diagonal  2nd:  P 
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Higher  Order  Moment  Conservation 


0 


9  Moments  defined  as  Integrals  of  VDF:  Q  =  f  Qnfdvi 
9  Discrete  Version:  nf  w^S(v^)  such  that  Q  =  Y1  w^Q/  Y2 
9  Merged  Particles  have  4  DOF  each:  w,  vx,  vz 
9  Number  of  Moments  Conserved  from  Number  of  DOF 


3  Particles  -  Mass,  Momentum,  Full  2nd:  P  & 
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Cartesian  Moments 


VxVj 

VxVy  iyty  Vy^ 
Wi  VyW 


Moment 

Order 

Mass 

Mass  Flux 

Momentum  Flux 

0th 

Ist 

2nd 

=  w 

=w  -Vi 

EwWvf'vf  |  W'W] 
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Higher  Order  Moment  Conservation 


0 


©  Moments  defined  as  Integrals  of  VDF:  Q  =  f  Qnfdvi 
©  Discrete  Version:  nf  w^5(v^)  such  that  Q  =  w^Q/ 
©  Merged  Particles  have  4  DOF  each:  w,  vx,  vy,  vz 
©  Number  of  Moments  Conserved  from  Number  of  DOF 


Cartesian  Moments 


Moment 

Order 

Mass 

0th 

^2  w&>  =  w 

Mass  Flux 

Ist 

=  w  ■  vj 

Momentum  Flux 

2nd 

=w-Wj 

Energy  Flux 

3rd 

Y.wtovPivW)2  = 

4  Particles  -  Mass,  Momentum,  Full  2nd ,  Energy  Flux:  qi 
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Extension  to  Hybrid  Methods 


Merge  Quantities  Needed  for  Hybridization: 

®  Reconstructed  VDF  Natural  extension  to 
Fokker-Planck/Boltzmann  Solvers 

o  Higher  Moment  Merges  would  Facilitate 
extension  to  Hybrid  Euler,  Navier-Stokes, 
13 -moment,  and  Beyond 

o  Reversal  of  VDF/Moments  to  Particles 
would  Enable  Particle  Generation  in 
Transition  Zones 


□  ►  <3  gp  >  <  ~  ►  <  -5  >  ^  '0<\0' 


Robert  Martin  (AFRL/RZSS) 


Distribution  A:  Approved  for  public  release;  unlimited  distribution 


12/13 


End 

Thank  You 

Questions? 
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